%%  mf_analyzer for Storage Tank for the visuals
% todo:
% - movie, lots of platform dependencies for now only option 3 and 4 work
% - plot streamlines

% b.f.destombe@student.tudelft.nl

close all; clear all; clc;

n = 1;    % setup which number you want to visualize
load setup setup
if isfield(setup(n),'outputPath'), saveOutput=true; end    % outputpath is in '.\output\current time'
%%  Options
% remember to set saveOutput in Runall if you want files to be saved
plotWellConc    = 0;    % plot concentration well and quantities
plotSetup       = 1;    % plots the aquarium, needed if you plot the 3D concentrations or streamlines
plotConc        = 1;    % plots the concentrations
plotStrmLines   = 1;    % plots streamlines
plotIsoPsi      = 0;    % plots surfaces with same pressureheads
checkBudget     = 0;    % Debug, check quantities
checkHeads      = 0;
makeMovie       = 0;    % Stub, 0=off;1=normal[not working], renders on the screen aswel;2=fast, renders faster and generates no plots, just the avi.[not working]; 3=make png figure strip; 4=eps figure strip; remember to set saveOutput in Runall
dim2D           = 1;    % front view
equiConc        = 1;
stub            = 0;

%% preperation, read datafiles
if (plotWellConc || plotConc || plotStrmLines || equiConc)
    C   = readMT3D('MT3D001.UCN');
    C   = maskHC(C,[0 1],[0 1]);
    for i   = 1:length(C)   % convert concentrationvalues to added-salt-content [kg/m3]
        C(i).values     = C(i).values.*(1075-1004);
    end
end

if (plotStrmLines || plotIsoPsi)
    B_  = readBud([setup(n).basename,'.BGT']);B=mf_Psi(B_);
    
    %     maskfunction to mask within boundaries; maskBud masks outside of
    %     boundaries.
    B   = maskBud(B,50);B=maskBud(B,-50);
    threshold= 0.999;  % to prevent roundoff errors
    MASK     = 0.001;
    for it=1:length(B)
        for iLbl=1:length(B(it).label)
            B(it).term{iLbl}(abs(B(it).term{iLbl})<threshold*MASK)=0;
        end
    end
end

if checkHeads
    H   = readDat([setup(n).basename,'.hds']);
end

%%  plotWellConc
if plotWellConc
    for i = 1:length(setup(n).well)     %not tested for multiple well structures
        setup(n).well    = setWellCout(setup(n).well,C,i);
        f = figure();
        hold on
        subplot(2,1,1); setup(n).well.plotCout(i);   title('Concentration [ ]')
        subplot(2,1,2); setup(n).well.plotQ;         title('Quantities [cm3/min]')
        drawnow;
        if saveOutput
            print(f,'-depsc',strcat(setup(n).outputPath,setup(n).basename,'plotWell_',num2str(i),'Series_',num2str(setup(n).setupVars.SERIES)));
        end
    end
end

%%  check Budget
if checkBudget
    zonebudget(Braw);
end

%%  check heads
if checkHeads
    Hcell   = struct2cell(H);
    Hmat    = cell2mat(Hcell(8,:));
    fprintf('The minimum head value is %f [cm]\n', min(Hmat(:)));
    fprintf('The maximum head value is %f [cm]\n', max(Hmat(:)));
end

%%  Plot Setup in 3D
if plotSetup
    %     Create grid
    if makeMovie == 2,	fig = figure('visible','off');  hold on;  % create handle
    else                fig = figure('visible','on');   hold on;  end
    
    Ax1=subplot(1,1,1);
    %     setup(n).gr.plotMesh;
    
    if makeMovie ~= 4,  lighting phong; end     % The rendermachine that can output eps file, does not support lighting
    
    %     plot aquarium
    glass.thickness = 0.5;	glass.alpha = 0.05;	glass.color = [0.1 0.1 0.8];    % Glass properties
    plotcube([-max(setup(n).gr.xGr) -max(setup(n).gr.yGr) -glass.thickness],      [max(setup(n).gr.xGr) max(setup(n).gr.yGr) min(setup(n).gr.zGr)], glass.alpha, glass.color);    % Bottom
    plotcube([-max(setup(n).gr.xGr) glass.thickness abs(min(setup(n).gr.zGr))],   [max(setup(n).gr.xGr) max(setup(n).gr.yGr) min(setup(n).gr.zGr)], glass.alpha, glass.color); % Back
    plotcube([glass.thickness -max(setup(n).gr.yGr) abs(min(setup(n).gr.zGr))],   [max(setup(n).gr.xGr) max(setup(n).gr.yGr) min(setup(n).gr.zGr)], glass.alpha, glass.color); % RHS
    %     injectionbarrier
    plotcube([-diff(setup(n).grExtra.xhInjB) -diff(setup(n).grExtra.yhInjB) abs(diff(setup(n).grExtra.zhInjB))], [max(setup(n).grExtra.xhInjB) max(setup(n).grExtra.yhInjB) min(setup(n).grExtra.zhInjB)], glass.alpha, glass.color); % horizontal
    plotcube([-diff(setup(n).grExtra.xvInjB) -diff(setup(n).grExtra.yvInjB) abs(diff(setup(n).grExtra.zvInjB))], [max(setup(n).grExtra.xvInjB) max(setup(n).grExtra.yvInjB) min(setup(n).grExtra.zvInjB)], glass.alpha, glass.color); % vertical
    %     plot well
    setup(n).well.plot3D('LineWidth',5);
    %     plot RECH; Not being used
    if 0
        RECHactive  = sum(abs(RECH),3);
        [iyRECH ixRECH izRECH]  = ind2sub(size(IBOUND),find(0 < RECHactive | RECHactive < 0));
        
        RECHplot    = scatter3(setup(n).gr.xm(ixRECH),setup(n).gr.ym(iyRECH),setup(n).gr.zm(izRECH));
    end
    
    %     plot HFB
    ixHFB   = setup(n).HFB(1,5);
    iyHFB   = setup(n).gr.yGr;
    izHFB   = 1:(max(setup(n).HFB(:,1))+1);
    
    p1  = [setup(n).gr.xGr(ixHFB) min(setup(n).gr.yGr) max(setup(n).gr.zGr(izHFB))];
    p2  = [setup(n).gr.xGr(ixHFB) max(setup(n).gr.yGr) max(setup(n).gr.zGr(izHFB))];
    p3  = [setup(n).gr.xGr(ixHFB) max(setup(n).gr.yGr) min(setup(n).gr.zGr(izHFB))];
    p4  = [setup(n).gr.xGr(ixHFB) min(setup(n).gr.yGr) min(setup(n).gr.zGr(izHFB))];
    x   = [p1(1) p2(1) p3(1) p4(1)];    y   = [p1(2) p2(2) p3(2) p4(2)];    z   = [p1(3) p2(3) p3(3) p4(3)];
    
    HFBplot     = fill3(x, y, z, 1);    alpha(HFBplot,0.5);
    
    %  plot properties
    xlabel('x [cm]');                 % set x-label
    ylabel('y [cm]');                 % set y-label
    zlabel('z [cm]');                 % set z-label
    c=colorbar('horiz');
    ylabel(c,'NaCl kg/m3')
    axis vis3d equal tight;
    
    if makeMovie < 2, drawnow;  end
    if dim2D, view([0 -1 0]);   end
end

if (3 > makeMovie > 0)
    aviobj=VideoWriter([setup(n).outputPath setup(n).basename '_setup_' num2str(n)]);   % create different video for each serie
    aviobj.FrameRate = 10;
    aviobj.Quality = 90; %1=Lowest quality, 100=Max quality
    open(aviobj);
end

%%  plot the action Concentration and Streamline
% create needed mesh
[U.X,U.Y,U.Z]=meshgrid(setup(n).gr.xc,setup(n).gr.yc,setup(n).gr.zc);

for i=1:setup(n).perdata.nper
    
    if (plotConc ~= 1 && plotStrmLines ~= 1 && plotIsoPsi ~= 1 && equiConc ~= 1),   break;  end
    
    hold on
    
    if plotConc
        if  i == 1
            [sliceH]    = slice(U.X,U.Y,U.Z,C(i).values,[setup(n).gr.xc(1) setup(n).gr.xc(end)],[setup(n).gr.yc(1)],[],'cubic');
            set(sliceH,'FaceAlpha',0.15,'FaceColor','interp','FaceLighting','gouraud','EdgeColor','none','DiffuseStrength',1)
        else
            cdata = sliceCData(U.X,U.Y,U.Z,C(i).values,[setup(n).gr.xc(1) setup(n).gr.xc(end)],[setup(n).gr.yc(end)],[],'cubic');
            for iii = 1:length(cdata)
                set(sliceH(iii),'CData',cdata{iii})
            end
        end
%         set(gca,'cLim',[0 72])
    end
    
    if equiConc
        if i == 1
            isorangeEquiConc=unique([linspace(0,72,8),0]);            
        end
        
        %	frontal view, over the y ax; will display flow from left to
        %	right
        for ii=1:length(isorangeEquiConc)
            [f,v] = isosurface(U.X,U.Y,U.Z,C(i).values,isorangeEquiConc(ii));
            hpatchEquiConc(ii) = patch('Faces',f,'Vertices',v,'FaceAlpha',0.35,'FaceColor','interp','EdgeColor','none');
            isocolors(U.X,U.Y,U.Z,C(i).values,hpatchEquiConc(ii));
        end
    end
    
    if plotStrmLines
        
        if i == 1
            isorangeStrm=unique([linspace(-8,8,20),-.1:.05:.1]);
        end
        
        V.X     = cumsum(B(i).term{2},3);
        V.Y     = cumsum(B(i).term{3},3);
        V.Z     = cumsum(-B(i).term{4},3);
        
        %	frontal view, over the y ax; will display flow from left to
        %	right
        for ii=1:length(isorangeStrm)
            [f,v] = isosurface(U.X,U.Y,U.Z,V.X,isorangeStrm(ii));
            hpatchStrm(ii) = patch('Faces',f,'Vertices',v,'FaceAlpha',0.35,'FaceColor','interp','EdgeColor','none');
            isocolors(U.X,U.Y,U.Z,C(i).values,hpatchStrm(ii));
        end
        
        [sliceStrmH averts] = streamslice(U.X,U.Y,U.Z,V.X,V.Y,V.Z,[],[setup(n).gr.ym(end)],[]); % Generate linedata for streamlines and
        sl2     = streamline(averts); % plot arrows
        set(gca,'cLim',[0 72])
        
    end
    
    fprintf('time plotted is %f\n',setup(n).perdata.time(i));
    
    set(gca,'DrawMode','fast')
    drawnow;
    
    for ii=1:length(isorangeStrm), delete(hpatchStrm(ii)); end
    for ii=1:length(isorangeEquiConc), delete(hpatchEquiConc(ii)); end

    
    delete(findobj(gca,'Type','line'));
    
    if (makeMovie == 1 || makeMovie == 2)
        h=getframe;
        
        if makeMovie == 1,  pause(.2);  end
        writeVideo(aviobj,h);
    else
        if (makeMovie == 3 && saveOutput), print(fig,'-dtiff','-opengl','-r150',strcat(setup(n).outputPath,'plot3D_time_',num2str(j))); end   % vectorplotting/eps cant handle a set alpha that wel..
        if (makeMovie == 4 && saveOutput), print(fig,'-depsc','-painters',strcat(setup(n).outputPath,'plot3D_time_',num2str(j))); end   % vectorplotting/eps cant handle a set alpha that wel..
    end
end

if (3 > makeMovie > 0), close(aviobj); end


if stub
    
    %     fprintf('injectionrate is %f\n',setup(n).well.Q(i));
    %     fprintf('extraction concentration is %f\n',setup(n).well.Cout(i));
    
    
    % for j = 1:max(setup(n).perdata.time) % slicecontour loop
    % for i = 1:length(setup(n).perdata.time);
    %     convert time to nper; find nearest time
    %     if j<min(setup(n).perdata.time), i=1; else i     = interp1(setup(n).perdata.time,1:length(setup(n).perdata.time),j,'nearest'); end
    
    
    [sx, sy, sz] = meshgrid([min(setup(n).gr.xm(:)),max(setup(n).gr.xm(:))], setup(n).gr.ym, setup(n).gr.zm);
    plot3(sx(:), sy(:), sz(:), 'bo', 'MarkerFaceColor', 'b');
    streamline(stream3(U.X,U.Y,U.Z,V.X,V.Y,V.Z,sx,sy,sz))
    
    [sx, sy, sz] = meshgrid([min(setup(n).gr.xm(:)), max(setup(n).gr.xm(:))], setup(n).gr.ym, min(setup(n).gr.zm(:)));
    plot3(sx(:), sy(:), sz(:), 'bo', 'MarkerFaceColor', 'b');
    streamline(stream3(U.X, U.Y, U.Z, V.X, V.Y, V.Z, sx, sy, sz))
end